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1. Introduction 

In the present contribution we shall analyze the effects of imbalancing the populations 
of two-component trapped fermions in the BEC (strong-coupling) limit of the attractive 
interaction between fermions of different components. In particular, we shall derive a set 
of coupled equations which describe composite bosons and excess fermions in this limit, 
starting from the gap equation with two different fermionic chemical potentials. Care will 
be used to include in these equations the processes leading to the correct dimcr-dimcr 
and dimer-fermion scattering lengths, which require us to consider beyond-mean-field 
effects. Numerical results will be presented for the density profiles of composite bosons 
and excess fermions, which are relevant to the recent experiments with trapped Fermi 
atoms. Results for the formation of vortex patterns in the presence of density imbalance 
will also be presented. 

The interest in imbalanced populations of fermions with different spins originated long 
time ago from condensed-matter physics, where a magnetic field can alter the populations 
of spin | and [ electrons [TJ [2] . 

Quite generally, density imbalance Sn — — ny between "spin up" and "spin down" 
fermions introduces a new degree of freedom in the system, and may thus possibly lead 
to the removal degeneracies and the occurrence of novel phases [3J H] . 

In this context, trapped cold Fermi atoms offer a unique possibility for observing 
the consequences of density imbalance. These systems has lately be object of mounting 
interest, both experimentally and theoretically, as they allow one to explore the BCS-BEC 
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crossover by controlling the interaction between fcrmions of different components via the 
use of Fano-Feshbach resonances. Recently, two experimental studies with imbalanccd 
populations have raised novel interest in these systems [6] . Density profiles of the two 
fcrmionic species as well as vortices have been detected. A quantum phase transition 
to the normal state on the BCS side of the crossover as well as phase separation in the 
crossover region have been identified. 

From a theoretical point of view, the BCS-BEC crossover gets modified as the many- 
body problem becomes richer in the presence of density imbalance. The interest in this 
problem has then involved not only cold-atom and condensed-matter physics, but also 
nuclear and subnuclear physics. On the BEC side of the crossover, it further gives one the 
opportunity of embedding into the diagrammatic structure the processes [7J [5] leading 
to the correct values of the scattering lengths for two composite bosons (as) and for a 
composite boson with an excess fermion (obf)- 

The effects of density imbalance on fermionic superfluids were originally studied in 
the weak-coupling (BCS) limit of the crossover both for the homogeneous [9] and trapped 
case [TO]. Only recently these calculations have been extended to cover the BCS-BEC 
crossover [TJJ , and to consider the effects of the trap [TU [T3] [TJ] which are essential to 
account for the experimental results with density imbalance. We thus begin by giving 
in the next Section a brief account of the mean-field approach for the homogenous case 
when different spin populations are considered, with emphasis to the strong-coupling 
(BEC) limit to which we shall eventually restrict in the later Sections when considering 
the trapped case. 



2. Mean- field treatment for the homogeneous case 

The microscopic BCS theory of superconductivity is most conveniently formulated in 
terms of fermionic single-particle Green's functions [15j . Due to the presence of sponta- 
neous broken symmetry, anomalous averages need be considered together with normal 
averages, leading to a 2 x 2 matrix for the Green's functions. 

This formulation can be readily extended to include different populations for the two 
fermionic species (labeled by spin-"f and spin- J,) which mutually interact via a contact 
potential. This is done by considering two different chemical potentials /if and fi^ for 
the two species, so that the equation of motion for the fermionic single-particle Green's 
functions Gy(k, uj s ) within mean field reads: 

f " 

1 



(1) 



-A* 



k 

2^7 -Vl 



Gn G12 
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Here, lo s = (2s+1)tt / (ksT) (s integer) is a fermionic Matsubara frequency at temperature 
T, k a wave vector, m the fermion mass, and A the gap function. The novelty introduced 
in eq. (JlJ) by population imbalance is the presence of two different chemical potentials in 
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the diagonal matrix elements on its left-hand side. This seemingly minor difference will, 
however, yield significant different results upon inverting eq. ([I]) in favor of G,j . 
From this inversion one obtains, in particular: 



n = n] +7H =J2 [e^GuQe) - **** 'Gzt(kj\ 
k 

dk f g(k) 
(27r)3 1 E(k) 



(2) = / 7^3 \ 1 " ^ [1 " f(E + m - f(E-(k))] 



for the total density (where r\ is a positive infinitesimal), 

<5n = n T - n; =2 [e^Gi^fe) + e~^' ) G 22 (A : )] 

/r/k 
^3{/(£+(k))-/(£?_(k))} 

for the density difference, and 

m _ /■ dk f l-/(£+(k))-/(£_(k)) m 



47ra F 7 (2tt) 3 [ 2£(k) k 2 

for the gap equation where the strength of the contact fermionic attraction has been 
replaced by the fermionic scattering length of by suitable regularization |16j . In these 
expressions, f(E) = (e^ E + l)" 1 is the Fermi distribution function with (3 = l/(k B T) 
and we have introduced the notation: 



(5) E ± (k) = E(k) ±8( , E(k) = ^(k) 2 + |A|2 
and 

(6) f(k) *_a+a , « = az£t 

Ira I I 

[only the case of equal fermion masses will be considered throughout] . Note from eq. 
that, in the low temperature limit, the only way to sustain a nonvanishing value of Sn is 
to have either E + (k) or E-(k) negative. For definiteness, we shall assume n-f > nj,. 

It is known from the theory of the BCS-BEC crossover that the dimensionless pa- 
rameter (hp a?) ~ 1 controls the evolution from the weak-coupling BCS regime (where 
of < and (fepap) -1 ~ —1) to the strong-coupling BEC regime (where ap > and 
(fep Of) _1 ~ +1), with the "crossover" region being limited in practice to the interval 
— 1 ~ (fcp of) -1 ~ +1. Here, kp is the Fermi wave vector related in the standard way to 
the total density. 

In the following, we shall mostly be interested in the strong-coupling (BEC) regime at 
low temperature. In this limit, one expects the presence of density imbalance to produce 
a density of composite bosons formed by pairing a fermion of spin f with a fermion 
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of spin I, plus a density Sn = rij — rij of excess fermions of spin |. Correspondingly, one 
finds from eqs. 

(7) /i| = - e - « T + fj, B 

where eo = rra 1 a 2 is the binding energy of the associated two-body problem and 

(8) /i T = e F {Sn) + ( 27rflBF ] nQ 



, n s _ / 47ra B \ / 27ra BF \ s 

(9) /i B = + on . 

\ m B J V m BF / 

In these expressions, e F (Sn) is the Fermi energy corresponding to Sn, n B plays the role of 
the chemical potential of the composite bosons, and no is the condensate density which 
equals (n — Sn)/2 under the present circumstances. In addition, m B = 2m is the mass of 
a composite boson, m BF = |m is the reduced mass of the two-body system made up of 
a composite boson and an excess fermion, while as = 2of is the scattering length for the 
scattering of two composite bosons and a B F = (8/3)of for the scattering of a composite 
boson and an excess fermion. These values for a B and a B F are specific to the mean-field 
treatment and correspond to the Born approximation for the scattering processes. 

After these preliminary considerations of general relevance in the presence of imbal- 
anced fermion populations, we pass now to consider the trapped case which is specifically 
relevant to the experiments with cold atoms. 



3. Mean-field treatment for the trapped case 

At a formal level the mean-field treatment for the inhomogeneous case proceeds along 
similar lines as for the homogenous case discussed in the previous Section, the difference 
being that the equation of motion for the fermionic single-particle Green's functions 
Gij(r,r';u> s ) is now given by: 



(10) 



-A(r)* 



-A(r) 
iu) s + Hi(r) _ 



G21 G22 



S(r - r') 



1 
1 



Here, Tt a (r) = — S— + V(r) — [i a is the single-particle Hamiltonian in the presence 
of the trapping potential V^r). Equation (fTU|) is the Green's function version of the 
Bogoliubov-de Gennes equations [17j . which are often used to describe inhomogeneous 
superconductors. 

As in the previous Section, we are mostly interested in the strong-coupling BEC regime 
at low temperature. In the case of equal fermion populations (whereby /i| = /i^), it has 
been shown [TS] that the Gross-Pitaevskii equation for composite bosons can be derived 
in this limit from the Bogoliubov-de Gennes equations (fT0|) . This mapping enables one to 
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exploit the results obtained directly from the more manageable bosonic Gross-Pitaevskii 
equation and, when needed, to use them as benchmarks for the fcrmionic calculation 
in the limit. This analysis has recently been extended to the imbalanced case when 
MT A*l- I n this case, in the place of the Gross-Pitaevskii equation one ends up with 
two coupled equations that describe the simultaneous presence of composite bosons and 
excess fermions [T2"] . As the analysis proceeds along similar lines in both (balanced and 
imbalanced) cases, we indicate it here schematically for convenience. 

One starts by rewriting the Bogoliubov-de Gennes equations (flC)]) in the integral form: 



(11) G(r, r';uj s ) = &(r, r'; w.) + j dv" G (r, r"; u a ) B(r") G(r", r'; u.) 



So(r,r';w s |^) 

-^ (r',r;-w s |/i;) 



where 

(12) Go(t,v';uj s ) = 
is the matrix of the noninteracting Green's functions which satisfy the equation 

(13) [iuJ s - H a {T)\ g (r,r';uM = <5(r-r'), 

being subject to the same trapping potential V(r) entering the single-particle Hamilto- 
nian Ti a (r). In addition, the matrix 



(14) 



B(r) = 



A(r) 
A*(r) 



contains the effects of the interaction via the gap function. 

At this point, expansion of fti(r,r';w s ) and £22 ( r , r'; u> s ) up to order A 2 yields the 
expressions for the local densities ^t( r ) an d n i( r )> m the order, while expansion of 
012 (r, t';uj s ) up to order A 3 yields an equation for the local gap A(r). The gap equation 
can be cast in the form [i"2l: 



(15) 



2m B 



$(r) 



2]/(r) + ^!fe(r) 



m 



$(r) + l^|$(r)| 2 $(r)=/x B $(r) 
m B 



where again me = 2m, a B = 2ap, a B F = |of, I^b = m + fJ>i + £0i while <i>(r) 

T "g° F A(r) plays the role of the bosonic condensate wave function. 
By a similar token, the density of excess fermions reads: 



(16) Sn(r) = n T (r) - n;(r) = 



dk 



(2tt) 3 j V 2m 



37ra B p 



■;?) 



|$(r)| 2 - m 



Note that the two coupled equations (fT5]) and (JTH)) embody the mutual effects of the 
bosonic distribution |$(r)| 2 and the fermionic distribution 6n(r). Finally, the expression 
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for the density is: 

(17) n(r) = ri T (r) + n x (r) = 5n(r) + 2 |$(r)| 2 

which completes the set of three coupled equations for <I>(r), Sn(r), and n(r). 

The remaining problem is that the scattering between composite bosons (as embodied 
by as) and between a composite boson and an excess fermion (as embodied by obf) has 
been treated so far at the lowest order within the Born approximation, corresponding to 
the values ae = 2of and Obf = f^F obtained by our derivation. This points to the need 
of going beyond the mean-field treatment to include the full set of scattering processes 
for Ob and gibf- In diagrammatic language, this implies identifying additional fermionic 
diagrams containing diagrammatic sequences for cib and obf as sub-units. 

4. Exact equations in the dilute case 

The validity of the equations obtained in the previous Section can be extended by 
improving on the values of the scattering lengths and obf ■ From the exact solutions (in 
real-space representation) of the three- [19] and four-body [20] problems it is known that 
the correct values are czbf = 1.18of and as = 0.6ap, in the order. These values have 
also been determined diagrammatically (in wave-vector representation) in the limit of 
vanishing density, for obf in ref. [5] and for as in ref . [7J . To improve on the derivation of 
the equations obtained in the previous Section, one has thus to embed the diagrammatic 
sub-units, which identify as and obf in the limit of vanishing density, into the many-body 
structure at finite density. This can be achieved as follows. 

As a natural extension of what is done for determining the condensate density of point- 
like bosons [21], in the BEC limit of interest here the gap equation can be interpreted 
as the condition of vanishing "tadpole" insertions for composite bosons. At the mean- 
field level considered in the previous Section, the diagrams representing this condition 
are depicted in fig. 1(a), with the understanding that a composite-boson propagator 
with zero four-momentum can be inserted from the left while the gap A corresponds 
to a condensate line. These diagrams account for the dimer-fermion and dimcr-dimer 
scattering processes within the Born approximation. The additional tadpole diagrams, 
which yield the correct values of obf and as, are depicted in figs. 1(b) and 1(c), in 
the order. Here, the dimer-fermion and dimer-dimer scattering processes exclude the 
Born-approximation contributions already included at the mean- field level in fig. 1(a). 

With this procedure, one ends up formally with the same coupled equations (fl5|) - 
(|17|) which, however, now contain the correct values of as and obf- In this context, a 
comment is worth on how to factor out from the diagrams of fig. 1(b) the product obf 
times 5n that enters eq. (| 1 5[) . The point is that the integration over the wave vector P 
is bounded within the Fermi sphere of radius */2m^, while the remaining integrations 
over q, q', • • •, extend outside this Fermi sphere. One may, accordingly, neglect the 
P-dependence everywhere in the diagrams of fig. 1(b) except in the fermion propagator 
labeled by P and corresponding to spin-| fermions. The density of excess fermions results 
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Fig. 1. - Tadpole diagrams for composite bosons describing the gap equation in the BEC limit, 
(a) Mean-field contributions; Contributions beyond the Born approximation which include the 
scattering processes yielding the correct values of the (b) dimer-fermion and (c) dimer-dimer 
scattering lengths. Full lines represent the fermion propagators of given spin and broken lines 
represent the bare fermionic interaction. [Reproduced from ref. |12j.] 



in this way from the P-integration, with the remaining parts of the diagrams yielding 
the exact dimcr- fermion scattering matrix T 3 (0). This excludes, by definition, the Born 
contribution resulting from mean field. A similar analysis can be carried out for the 
density equation (| 16[) , where only the value of the dimer-fermion scattering length obf 
requires corrections beyond mean field [P^] . 

5. Numerical results and comparison with experiments 

We pass now to report on the numerical solution of the coupled equations (| 1 5[) - (flT|) 
with the exact values of ob and obf- By their very derivation, these equations are 
expected to be valid in the strong-coupling (BEC) regime where of > 0, which however 
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Fig. 2. - Density profiles 5n(r) (full lines) and no(r) (broken lines) vs. r for (/ifof) -1 = 3 and 
a — 0.2 (upper panel), a = 0.5 (middle panel), a = 0.8 (lower panel). The insets show the 
results when (fcpOF) -1 = 1 for the same values of a. Here, r is in units of Rtf and densities are 
in units of (iVf + Ni)/R TF . The Thomas- Fermi radius Rtf and the Fermi wave vector k F are 
defined for equal populations. [Reproduced from ref. |12].] 



extends as down as to (&f of) -1 ~ +1 for all practical purposes. 

We determine the density profiles Sn(r) and no(r) = |$(r)| 2 in the Thomas- Fermi 
(LDA) approximation (which corresponds to neglecting the kinetic energy term in eq. fl5])) 
for a spherical trap, as functions of the asymmetry parameter a — (JV-j- — A/j)/(iVf + N±) 
(where < a < 1). Numerical results [12] are shown in: 

(i) fig. 2 for the density profiles Sn(r) and n (r) vs. the distance r — |r| from the 
center of the trap, for three characteristic values of a and for the coupling (fcpap) _1 = 3 
on the BEC side of the crossover. The insets show the results for the smaller coupling 
(k-pa-p)^ 1 = 1. A spatial separation results between the condensed composite bosons 
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Fig. 3. - Column density profiles 5n(p) (full lines) and no(p) (broken lines) vs. p for (fear)' 1 = 3 
and a = 0.2 (upper panel), a — 0.5 (middle panel), a = 0.8 (lower panel). The insets show the 
results when (fcaF) -1 = 1 for the same values of a. Here, p is in units of -Rtf and densities 
are in units of (iV T + iVj ) / i?^F • 



and the excess fermions, which appears sharper for the smaller coupling, correspond- 
ing to enhanced effects of the dimer-fermion repulsion. Upon approaching unitarity 
[(fcpap) -1 = 0], there thus appears a tendency toward phase separation with a super- 
conducting core of fully paired fermions surrounded by a cloud of excess fermions. For 
the couplings here considered, this phase separation occurs for all values of a. Note 
that, for each value of a, the maximum of 5n(r) occurs where n (r) vanishes, and that 
there occurs a progressive size shrinking of no(r) for increasing a, with a simultaneous 
penetration of 5n(r) toward the center of the trap. 

(ii) fig. 3 for the corresponding column density profiles n(p) — Jdzn(p, z) vs. p = 
\J x 2 + y 2 . Phase separation appears now less visible even for [kpa-p)^ 1 = 1, since n(p) 
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0.2 0.4 0.6 0.8 

r/Rp 

Fig. 4. - Column density profiles 5n(p) (full lines) and no(p) (broken lines) vs. p when (fcpap) -1 = 
1 and a = 0.5, as obtained with the exact (E) and mean-field (MF) values of the scattering 
lengths. Units are the same as in fig. 3. 



leaks toward p = where it acquires a finite value. 

(iii) fig. 4 for the comparison of the column density profiles when (fcpap) -1 = f and 
a = 0.4, as obtained with the exact and mean-field values of the scattering lengths. The 
use of the exact values of the scattering lengths results in density profiles which are more 
compressed toward p = 0. This feature can be conveniently characterized by the values 
of the critical radius R c where the condensate vanishes. 

(iv) fig. 5 for the critical radius R c of the condensate density, which is plotted vs. a 
for different couplings and using the exact (a) or mean-field (b) values of the scattering 
lengths. This quantity identifies also the position of the maximum of the density of excess 




0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 



a a 

Fig. 5. - Critical radius R c of the condensate density vs. the asymmetry parameter a (normalized 
to the value -R c (0) at a = 0) for three coupling values, obtained by using the exact (a) or mean- 
field (b) values of the scattering lengths. [Adapted from ref. [12].] 
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0.8 1.2 1.6 

a BF /a F 

Fig. 6. - Ratio v between the column density of excess fermions at the center of the trap and 
its maximum value at R c vs. the ratio aBp/aF, for a = 0.5 and three coupling values. [Adapted 
from ref. [12].] 



fermions. Note the inverted sequence of the curves corresponding to the different values 
of (fcpaF)" 1 , as calculated with the exact or with the mean-field values of the scattering 
lengths. This feature could be subject to experimental verification. 

(v) fig. 6 for the ratio v between the column density of the excess fermions at the 
center of the trap and its maximum value at the critical radius R c vs. obf/of for a = 0.5 
and different couplings. The use of column density to obtain the ratio v stems from our 
finding that column density profiles have a more marked dependence on obf- Indeed, 
the marked dependence on obf should make it possible to extract the expected value 
1.18 of obf/af from the experimental data, using the plots of fig. 6 as calibration curves. 
This procedure is indicated schematically in the figure. 

Note that, although the results presented in this Section have been obtained for an 
isotropic (spherical) trap, they can be also utilized for an anisotropic (ellipsoidal) trap, 
as indicated in Appendix A. 

Comparison with the experimental results by the MIT [5] and Rice [6] groups for the 
density profiles ot the two imbalanced fermion species can be established at this point, 
at least as far as the BEC side of the crossover region is concerned. In this regime and at 
low temperatures, the experiments indicate that the dominant effect of density imbalance 
is to phase separate the system in two components, with an inner superfluid region where 
fermions of different species balance each other and an outer normal region where the 
excess fermions reside. This effect was correctly reported for the first time in ref. [T^] 
and it is clearly evidenced in the above figures. 
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6. Extension to vortices (rotating frame) 

Unambiguous detection of the superfiuid phase for trapped Fermi atoms with balanced 
populations has been achieved by the observation of vortex lattices when the system is 
put into rotation [22] . By a similar token, in the presence of density imbalance between 
the two fermionic species, detection of the superfiuid region has relied on the observation 
of vortices in that region [5] . 

To deal with vortex patterns originating from fermions in the presence of density 
imbalance, we extend the previous treatment holding on the BEC side of the crossover 
in the following way. We begin by noting that the two coupled equations (fT5j) and (fTB|) 
for the condensate wave function <&(r) and the density of excess fermions Sn(r) could 
alternatively be derived from a suitable energy functional E[<f>, 5n]. When the system is 
further put into rotation, this energy functional reads: 

E[$,6n] = Z fdxdy (-^-|V$(r)| 2 + mui 2 r 2 |$(r)| 2 + — |$(r)| 4 
J [ 4m m 

«5/3 4/3 i 

- ft$(r)*L z $(r) + ^— <5n(r) 5 / 3 + -mtu 2 r 2 5n(r) 

20m 2 

(18) + ^2£|$(r)| 2 <5n(r) - -mQ 2 r 2 Sn{r) 

m 2 

Here, Q is the rotation frequency, L z the angular momentum operator of the composite 
bosons, r = (x,y 7 0), while for calculation convenience the original ellipsoidal trap has 
been replaced by a cylinder of height Z with a harmonic radial potential. Note that in 
eq. (|18p the excess fermions in the rotating frame have been treated within a semi-classical 
description. Let, in fact, 



(19) n R (r) 



dp 1 



be the semi-classical fermionic distribution in the rotating (R) frame, where e#(p,r) = 

2 

f— + U(r) — SI ■ 1 the associated semi-classical dispersion with angular momentum 1. Since 
the last expression can be manipulated as follows 

P 2 

ej j p,r =f- + U(r) -p-(nxr) 
(20) = ^ + U(t) - \m{n x r) 2 

where p' = p — m£l x r in the last line, eq. (JH]) can eventually be cast in the form: 
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X/Rp X/R F 

Fig. 7. - Pattern of an isolated vortex in the balanced situation (a — 0) at the lower critical 
frequency for coupling values (a) (kyaF)^ 1 = 1 and (b) (kpap) -1 = 4. 



which agrees with the expression for Sn(r) one obtains from eq. (|18|) . 

Sticking to this equation, one has the choice of either solving the associated coupled 
differential equations for $(r) and <5n(r) or minimizing directly E[<f>, Sn] by varying <I>(r) 
and <5n(r). In practice, a "mixed" approach may sometimes be preferred 23J. Note that 
in the present treatment the Thomas-Fermi (LDA) approximation is altogether avoided. 
[We have also verified that this refinement leads only to minor changes for the density 
profiles presented in Section 5 in the absence of rotation.] 

We present numerical results for the parameters of the MIT trap [5] for which N = 
2 x 10 6 , uj± — 57Hz, and uo z = 23Hz. We begin by showing in fig. 7 the pattern of an 
isolated vortex in the balanced (a = 0) situation at the lower critical frequency f2 Cl when 
the first vortex enters the system, for the two cases (fcpap) -1 = 1 and fl Cl = 3.95Hz 




Fig. 8. - Triangular-shape vortex lattice in the balanced situation (a = 0) when Q = lQHz and 
for coupling values (a) (fcpap) -1 = 1 and (b) {k^a-p) -1 = 4. 
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Fig. 9. - Triangular-shape vortex lattice in the imbalanced situation with a — 0.4, £1 — 19Hz, 
and coupling values (a) (fcpap) -1 = 1 and (b) (SifOf) -1 = 4. The upper (lower) plots refer to 
the density of condensed composite bosons (excess fermions). 



(a), and (k-pciF^ 1 = 4 and fl ci — 6.11Hz (b). When the rotation frequency increases 
further, this pattern evolves in a triangular-shape vortex lattice shown in fig. 8 when 
ft = 19Hz for the same coupling values as in fig. 7. Such triangular pattern persists 
in the imbalanced situation, as shown in fig. 9 for a = 0.4 and the same frequency 
= 19Hz, when (a) (fepap) -1 = 1 an d (b) (k-pap) _1 = 4. In fig. 9 the upper plots 
picture the density of condensed composite bosons while the lower plots refer to excess 
fermions. Note that the excess fermions tend to pile up in the vortex cores, taking thus 
advantage of the fact that the density of composite bosons is depressed there. Further 
study along these lines is in progress [25] , 

7. Perspectives and open problems 

In this contribution, we have presented a theoretical account of the effects produced 
by imbalancing the populations of trapped fermions. Special emphasis has been placed 
to the strong-coupling (BEC) limit, where theoretical studies can be extended beyond 
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BCS mean field by exploiting the diluteness condition of the system. Several effects were, 
however, left out from our analysis. 

First of all, it appears relevant to combine density imbalance with mass imbalance of 
the two fermionic species, which may act to promote still novel phenomena in the system 
and it is being subject to experimental investigations at present. 

The present approach was confined to the BEC region where (fepap) -1 ~ +1. As 
the experiments [5J [5] actually span the whole crossover region -1 ~ (fc F ap)" 1 ~ +1, 
theoretical calculations should cover this region, too. In this context, a mean-field ap- 
proach could readily produce a qualitative overview of the effects taking place across 
the crossover region, including both density and mass imbalance. As in the case of the 
BEC region, however, mean-field calculations may lead to quantitative incorrect results. 
For instance, at unitarity and for T ~ 0, mean-field calculations [Hj result in the value 
a\ lF ~ 1 for the critical asymmetry parameter a c past which the superfluid region dis- 
appears from the system, while the experiment [5] yields a^, xp ~ 0.70. This is a clear 
indication that inclusion of pairing fluctuations beyond mean field is relevant in the 
crossover region to account for the experimental data. 

At higher temperatures, density imbalance may give access to "precursor" pairing 
(pseudo-gap) effects for the excess fermions, which could thus be evidenced in the normal 
phase even below the critical temperature T c . This may result in a strong suppression of 
thermal fluctuations, leading to a possible detection of an underlying Quantum Critical 
Point. 

Finally, we mention that a compelling test for theories would result from calculating 
the dependence of T c on a, as the experiments were able to determine the temperature 
of the system directly from the "tail" of the density profiles of the excess fermions |24) . 

Appendix A. 

Appendix A: Mapping of the anisotropic onto the isotropic problem 

The calculations presented in Section 5 refer to an isotropic harmonic potential (spher- 
ical trap). In the experiments with cold Fermi atoms, however, anisotropic harmonic 
traps are most often used. The question thus naturally arises to what an extent the the- 
oretical calculations done for a spherical trap can account for the experimental situation. 

The answer to this question relies on the validity of the Thomas-Fermi (LDA) ap- 
proximation, which applies when the number N of trapped atoms is large enough and 
was explicitly used in the calculations of Section 5. [As a matter of fact, when N is large 
enough, avoiding this approximation leads only to minor changes of the density profiles, 
as mentioned in Section 6.] 

This is because, within LDA, physical quantities acquire an r— dependence only 
through the local chemical potential 

(A.l) Kx,y,z) =/i--m^(i 2 + X 2 y y 2 + \\z 2 ) . 

Here, the anisotropy of the harmonic potential is specified by the parameters X y and A z . 
In particular, the density profiles of interest can be expressed as n(x, y, z) — n[/i(x, y, z)\. 
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To establish the desired mapping with the isotropic case, we introduce an equivalent 
isotropic problem with a total number of atoms N- lso = X y X z N, the harmonic frequency 

1/3 

ujiso = and the same Fermi energy Ep° = Ep = (2>u) x uiyU) z N) 1/0 of the original 
anisotropic problem. To determine the corresponding chemical potential /ii SO for this 
problem, we argue that the ratio pi so / 'Ep must be a universal function of the dimension- 
less quantities (kpap)~ 1 and T/Ep (where, by definition, kp — y / 2mEp). This implies 
that the ratio p,i so /Ep can depend on N only via Ep. As a consequence, pi so of the 
equivalent isotropic problem must coincide with the chemical potential p of the original 
anisotropic problem. 

For the density profile n{x,y, z) of the original anisotropic problem, this equivalence 
eventually implies that it can be expressed in terms of the density profile n>i so (x, Y, Z) of 
the equivalent isotropic problem as follows: 

(A. 2) n(x,y,z) = n l: 

where Y = X y y and Z = X z z. 



1 o , o 2 
2 "' 



Appendix B. 

Appendix B: Axial density profiles 

The density profiles reported in Section 5 were obtained within the Thomas-Fermi 
(LDA) approximation. This approximation tends to sharpen the density profiles at their 
edges, but is otherwise appropriate when the number of trapped atoms is large enough. 
From those results we have also concluded that the dominant effect of density imbalance 
is to phase separate the system into two components, an inner superfluid region and an 
outer normal region with excess fermions only. 

In this context, an argument was proposed by De Silva and Mueller [25] to evidence 
the occurrence of phase separation in the density profiles when the trapping potential is 
harmonic. Their argument relies on the validity of the LDA approximation and applies 
to the axial density n\(z), obtained from the density through a double integration: 

//>oo />oo 
dx dy n(x,y, z) = tt / dp n(p + z ) = tt / d£n(£) 
Jo Jz 2 

where £ = p 2 + z 2 and n(£) > 0. 

From the last line of eq. (|B. 1[) one readily concludes that ua(z) is a decreasing function 
of z. In addition, when n(x,y, z) presents a "hole" in the core region (say, for r < ro), 
tia{z) is constant in that region. This situation is pictured in fig. 10, where schematic 
plots of the radial density n(Q with a hole in the core region (upper panel) and of the 
corresponding axial density nx{z) with a flat behavior in the core (lower panel) are 
shown. This situation is a fingerprint for the occurrence of phase separation and applies 
to the density profiles shown in the insets of fig. 2. 

* * * 

We thank Stefano Simonucci, who performed the numerical calculations on vortices 
and prepared figs. 7-9. This work was partially supported by the Italian MIUR (contract 
Cofin-2005 "Ultracold Fermi Gases and Optical Lattices" ) . 
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Fig. 10. - Schematic plot of the density n(Q with a hole in the core region and of the corre- 
sponding axial density ua(z) with a flat behavior in the core. 
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